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In this article, we outline the modeling of a realistic dynamical model for comprehensive 
description of high energy heavy ion collisions. Comparing theoretical calculations and exper- 
imental data at RHIC, we give detailed discussions on the key ingredients for the construction 
of a multi-module model: initial condition, hydrodynamical expansion, hadronization, and 
freezeout processes. 
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§1. Introduction 



Since the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Lab- 
oratory (BNL) started its operation in 2000, a lot of discovery has been made and a 
lot of insight related to quantum chromodynamics (QCD) phase transition and the 
Quark-Gluon Plasma (QGP) has been gained. One of the most physically interesting 
and surprising outcomes at RHIC is the production of the strongly interacting QGP 
(sQGP). This accomplishment was realized by combining investigations from both 
experimental and theoretical sides. Because the QGP had been believed to be a 
weakly interacting system like the ideal gas, the discovery of sQGP opened up a new 
\ paradigm for the understanding of high temperature and/or high density QCD. 

Q^ ' The highlights of the RHIC experiments were (i) strong elliptic flow, which sug- 

gests early thermalization and early formation of collectivity; (ii) strong jet quench- 
. \ ing, which confirms the formation of hot and dense matter in collisions; (iii) con- 

. of deconfined hot quark soupP^ From studies of these experimental results with 

relativistic hydrodynamical models, jet energy loss mechanism, and recombination 
models, the discovery of sQGP at RHICP 1 was accomplished. Furthermore, heavy 
ion collision operation at the Large Hadron Collider (LHC), whose collision energy 
is around 15 times as large as that at RHIC, started in 2010. Such a high energy 
collision experiment gives us an opportunity to perform further investigation on the 
QCD diagram. 

Figure [1] shows a schematic sketch of the time evolution in relativistic heavy 
ion collisions on the basis of present understanding. When two heavy ions are ac- 
celerated to high energy, an extreme state which is described by static color charge 
and classical gluon fields, called the color glass condensate (CGC) state, is realized 
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inside each heavy ion. At the collision of the two heavy ions, high energy prompt 
photons and Drell-Yan dileptons, and jets are produced by hard scattering of quarks 
and gluons. After the collision thermalization is achieved in a short time. During 
that time, prethermal photons and dileptons are created and the entropy of the fire- 
ball increases. If the collision energy is large enough, QGP is also produced at this 
stage. Then hydro dynamical expansion starts, and a lot of interesting physical pro- 
cesses occur: collective flow formation, jet quenching, production of thermal photons 
and dileptons, and so on. As the temperature and density of the fireball decreases, 
hadronization of the QGP phase takes place. To understand the hadronization mech- 
anism in relativistic heavy ion collisions, the recombination/coalescence mechanism 
and the fragmentation mechanism are employed. Gradually the mean free path be- 
tween hadrons becomes large, and eventually freezeout happens. At this stage, final 
state interactions play an important role. 




Initial fluctuation hydrodynamic model final state interactions 



Fig. 1. Schematic sketch of relativistic heavy ion collisions. The wavy lines stand for photons and 
the arrows stand for jets. 

The present understanding of heavy ion collisions strongly suggests (Fig. [T]) that 
a multi-module modeling is indispensable for the description of entire history of 
heavy ion collisions. Knowledge of dominant physics at each stage has been accu- 
mulated, but a comprehensive model is still missing. For the construction of such 
a multi-module model, hydrodynamical models are a promising starting point, be- 
cause at present it is considered as one of the most reliable and successful dynamical 
models for understanding experimental data at RHIC and LHC comprehensively, 
especially for the description of the QGP phase. At the same time, it is easy to 
implement the latest developments in the physics of heavy ion collisions such as fluc- 
tuated initial conditions, the lattice QCD inspired equation of state, recombination 
mechanism for hadronization, and final state interactions in freezeout processes into 
a hydrodynamical model, as we will describe later in detail. 

In Tabs. U and [II] current hydrodynamical models are listed. In line with the 
physical picture shown in Fig.[TJ in particular, we pick up and compare the following 
features of each hydrodynamical model: dimension of hydrodynamical expansion 
(dim), initial condition (IC), equation of state (EoS), and treatment of freezeout 
process. In addition, we also compare the numerical schemes in solving the rela- 
tivistic hydrodynamical equation and observables calculated in each model in the 
tables. The importance of utilizing proper numerical schemes in solving the rela- 
tivistic hydrodynamical equation will be discussed in Sec. 0] Recent development 
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in relativistic viscous hydrodynamical models is remarkable. From the point of view 
of multi-module modeling, however, the status of ideal hydrodynamical models is 
considered to be more mature. 

Table I. Ideal hydrodynamical models. In the table, we use the following abbreviation. IC: ini- 
tial condition, G: Glauber model, CGC: color glass condensate, MC-G: Monte Carlo Glauber 
model, MC-CGC: Monte Carlo CGC, 1QCD: lattice QCD inspired EoS, SPH: smoothed particle 
hydrodynamics, PPM: piecewise parabolic method, CE: continuous emission, Obs: calculated 
observables, and PD: particle distribution. 



Ref. 


dim 


IC 


EoS 


scheme 


freezeout 


Obs 


Ham<P 


3+1 


NeXus 


Bag model 


SPH 


CE 


PD, V2, HBT 


Hirano 1 ^' 


3+1 


G, CGC 


Bag model 


PPM") 


cascade(JAM) 


V2 


Nonak<P 


3+1 


G 


Bag model 


Lagrange 


cascade(UrQMD) 


PD, v 2 


HirancPUU* 


3+1 


MC-G, MC-CGC 


1QCD 


PPM**) 


cascade(JAM) 


V2 


Petersen^ 


3+1 


UrQMD 


hadron gas 


SHASTA 


cascade(UrQMD) 


PD 


HolopainerP^ 


2+1 


MC-G 


1QCD 


SHASTA 


resonance decay 


V2 



Table II. Viscous hydrodynamical models. In the table, we use the following abbreviation. CD: 
central difference, and KT: Kurganov-Tadmor (KT) scheme. 



Ref. 


dim 


IC 


EoS 


scheme 


freezeout 


Obs. 


Romatschke^ 


2+1 


G 


1QCD 


CD 


single Tf 


V2 


Dusling^ 


2+1 


G 


ideal gas 




viscous correction 


V2 


LuzunP^ 


2+1 


G, CGC 


1QCD 


CD 


resonance decay 


V2 


SchenkfP 1 


3+1 


MC-G 


1QCD 


KT 


viscous correction 


V2, V3 


SonjF* 


2+1 


MC-G, MC-CGC 


1QCD 


SHASTA 


cascade(UrQMD) 


V2 


Chaudhurl^^ 


2+1 


G 


Bag model 


SHASTA 


viscous correction 


V2 


BozelPf 


3+1 


G 


1QCD 




THERMINATOR2 


VI, V 2 , HBT 



In this article, we outline the modeling of a realistic dynamical model for the 
description of relativistic heavy ion collisions in line with the physical picture shown 
by the schematic sketch (Fig. [1]). In the discussion, we refer to and compare with 
experimental data at SPS, RHIC, and LHC. The article is organized as follows. In 
Sec. [21 we review the initial condition for hydrodynamical models from the conven- 
tional Glauber type one to the latest attempt to include event-by-event fluctuated 
initial conditions. In Sec. [3j we present the basic concepts and ingredients of both 
ideal and viscous hydrodynamics such as relativistic hydrodynamical equation, equa- 
tions of state, and transport coefficients. We will also discuss interplay among jets, 
medium, and hydrodynamical expansion. In Sec.0J we review the numerical schemes 
which are listed in Tabs . HI and HTl and show the result of the newly developed scheme 
by one of the authors and her collaborators for a relativistic viscous hydrodynamical 

**' It is now customary in numerical hydrodynamics not to call a hydrodynamical computer 
program a PPM (piecewise parabolic methodp* one unless it fulfills not only parabolic interpolation 
of variables but also sharpening of discontinuity profiles and flattening of post-shock oscillations.^ 
Hirano's cod^^^^ executes only the first of the above three ingredients.^ Attention needs to 
be payed on this difference when comparing with other PPM codes and estimating its capability to 
capture shocks. 
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model. In Sec. [5j we explain the recombination model and show its utility in un- 
derstanding hadron observables and the QCD phase diagram. In Sec. [6J we explain 
chemical and thermal freezeout processes and discuss effects of final state interactions 
by comparing theoretical calculation and experimental data. Section [7] is devoted to 
summary and conclusions. 

§2. Initial conditions 

The hydrodynamical equations of motion requires inputs of initial conditions for 
all their dynamical variables, which are then evolved forward in time. These initial 
conditions are outside of the framework of hydrodynamical models and have to be 
determined by other means. Physically, they are determined by the process during 
the initial collision of the nuclei and the succeeding stage that makes the system 
approach to equilibrium, which is eventually reached at a time tq. Note that To, 
in principle, can depend on the coordinate space rapidity r/, while in practice it is 
assumed to be independent of rj. The equilibration time is a parameter since the 
equilibration mechanisms are still under debate and the first principle determination 
of the initial conditions of the equilibrated plasma phase has not been achieved EU'GSJ 

Historically, parameterized initial conditions for entropy density (or alterna- 
tively the energy densities) and the net baryon density have been usedF^'®'®'® 1 
In the transverse plane these distributions have been mainly parameterized based on 
Glauber-type models of nuclear collisions. In the longitudinal direction initial distri- 
butions inspired by Bjorken's scaling solution are often used. Then a few parameters 
remain to be fixed additionally in the initial condition, such as the maximum values 
of the energy or entropy density, and net baryon density. They are usually fixed 
by comparison with experimental data on single particle rapidity distributions and 
transverse momentum spectra. 

As a first trial, one can choose to set the initial longitudinal flow to Bjorken's 
scaling solutionis} and one can set the initial transverse flow to zero. This simplest 
initial flow profile has served as the basis for all further investigation. The possibil- 
ity of the existence of an initial transverse flow at tq was discussed, e.g., by Kolb 
and Rapp!^ The initial flow improves the results for P^-spectra an d reduces the 
anisotropy. This suggests that HBT analyses may be a sensitive tool for the deter- 
mination of the initial longitudinal flow. We note that, for the analysis of the final 
state longitudinal flow, the Yano-Koonin parametrization is effective. Hydrodynam- 
ical calculations during the early RHIC years did show serious disagreement with 
experimental data, especially for the ratio of R ou t/ Rside, leading to the notion of the 
HBT puzzle.® It turned out that the solution of the HBT puzzle is not so simple 
because it is related to all stages of hydrodynamical model: initial conditions, the 
equation of state, viscosity effect, and final state interactions. PrattP^ 1 showed that 
this HBT puzzle comes from not a single shortcoming of hydrodynamical models, 
but the combination of several effects; it is solved by mainly prethermal acceleration, 
a stiffer equation of state, and viscosity effect. 

Let us come back to the apparent early thermalization times found at RHIC. 
Usually it is argued that small initial times tq are needed to describe the elliptic 
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flow data as the elliptic flow builds up at the earliest stage of expansion when the 
eccentricity of the fireball is largest J3SJJ3D However, we note that with suitable sets 
of initial conditions and freezeout temperatures in fact a larger initial proper time 
is also compatible with data. Luzum and Romatschke show that three very differ- 
ent sets of the initial and freezeout temperature (Tj,Tf) — (0.29,0.14) GeV with 
t = 2 fm, (0.36,0.15) GeV with r = 1 fm, and (0.43,0.16) GeV with r = 0.5 fm 
- provide almost identical differential elliptic flows in their viscous hydrodynam- 
ical calculation.® This suggests that better constraints on initial conditions are 
indispensable to avoid incorrect conclusions from comparisons of hydrodynamical 
calculations with experimental data. 

Other approaches have been also taken in generating initial conditions. Color 
glass condensate-inspired initial conditions are becoming increasingly popular (Tabs. 
UlandHH). They feature larger eccentricities of the initial energy profile than Glauber- 
based models, which has significant implications on elliptic flow.® In these models 
additional dissipation during the early quark-gluon plasma stage is needed in order 
to achieve agreement with experiments.®'® Others models include the string rope 
model® and the pQCD + saturation model.'® In the latter the initial time tq is 
given by the inverse of the saturation scale, which is very small, i.e., r =0.18 (0.10) 
fm at RHIC (LHC). 

More recently there is a push to implement effects of event-by-event fluctuations 
in the initial conditions. In the NEXSPHERIO hydro model each event is created 
by the event generator NeXusP First they found that the existence of fluctuation 
in initial conditions improves the behavior of the elliptic flow as a function of the 
rapidity in hydrodynamical calculation.® They showed that the two artificial bumps 
in the elliptic flow as a function of the rapidity^* disappear if they take into account 
of initial fluctuations. 

An interesting observable, "Mach-Cone-like structure" , which is an angular cor- 
relation with the leading jet particle was reported at RHICP^'^ At first the origin 
of the structure is considered as the remnant of the interactions between jets and 
medium ([18]-[26] in RefP*). Models succeeded in giving the qualitative interpreta- 
tion of the structure, but they are not so successful in quantitatively describing the 
experimental data. For the first time in Ref!^* event-by-event fluctuations in the 
initial state were considered as a possible origin of the structure. A breakthrough 
of clear understanding of the Mach-Cone-like structure was done by detailed exper- 
imental analyses of triangular flow and higher harmonics. Current understanding 
is that the Mach-Cone-like structure is dominated by the triangular flow (See for 
example, RefP^). The triangular flow and higher harmonics are the coefficients in 
the Fourier expansion of particle yields as a function of the azimuthal angle (ft, 



— — oc l+2?;icos(0-6'i)+2t;2Cos2(0-02)+2'y3Cos3((/>-03)+2?;4Cos4((/>-6'4)+--- , 
dyacft 

(2-1) 

where v±, t>2, i>3, and v± are directed, elliptic, triangular, and quadrangular flows, 
respectively. 0i(i = 1,2, ... ) are constants and in principle independent from each 
other. The origin of the triangular flow and higher harmonics is fluctuation in initial 
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conditions 401 (Fig. [2|). In particular, v% vanishes, if the system starts with a smooth 
almond-shaped initial stated (Fig. [3]). Investigation of the relation between initial 
geometry and higher harmonics has been also carried out! 41 ^ For more quantitative 
analyses, however, contributions of final state interactions should be evaluated. 





Fig. 2. Elliptic and triangular flows from fluctu- 
ated initial state. 



Fig. 3. Elliptic flow from a smooth almond- 
shaped initial state. 



Here we make a comment on event-by-event fluctuations in the initial conditions 
in hydro dynamical models from the point of view of the numerical solution of the 
hydrodynamical equation. When the hydro dynamical simulation is performed with 
initial conditions with the event-by-event fluctuation, shock-wave capturing schemes 
should be used to describe the hydrodynamical expansion correctly. Otherwise, 
the effect of the fluctuation is smeared. It is known that most of the numerical 
schemes used in the calculations of the time evolution of the quark-gluon plasma do 
not satisfy this requirement. In other words, they introduce numerical viscosity to 
a non- negligible extent. In the early stage of the RHIC operation, hydrodynamical 
models gave us certain evidence of the existence of strongly interacting QGP at RHIC 
and are eventually regarded as the most reliable dynamical models for high energy 
heavy ion collisions. However, recent high statistical experimental data impose more 
rigorous numerical treatment on the hydrodynamical models. We will discuss this 
issue in detail later in Sec. HI 

2.1. Experimental Data and Discussion 

Among various kinds of experimental observables, thermal photons are one of the 
most promising ones for the investigation of the initial conditions of hydrodynamic 
models. PHENIX Collaboration reported that the excess from the superposition 
of pp collisions in its direct photon measurement at low transverse momentum Pj>. 
This excess follows the exponential distribution, which suggests that the thermal 
equilibrium is achieved! 4 ^* They extract the inverse-slope parameter from the photon 
spectra in Au+Au ^fs = 200 GeV, T = 221 ± 19 stat ± 19 syst . This value is larger than 
the critical temperature of the phase transition (~ 160 MeV), but it is smaller than 
the initial temperatures of hydrodynamic models at RHIC, T ~ 300 - 600 MeV. 
The value is interpreted as the average temperature in the whole hydrodynamic 
expansion. Contributions from the hadron phase are also important. 

It will be appropriate to point out here the importance of the hydro+micro 
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model as an instrument that relates the initial state condition to the final state 
observables. Very often, for example, final state multiplicities are directly predicted 
from the initial condition from the color glass condensate picture. For example, 
when the first LHC multiplicity data came outp" deviation of the experimental 
data from the prediction of the color glass condensate picture was found and was 
regarded as a serious problem. In fact, the system goes through several processes in 
which entropy and multiplicity change, such as the thermalization of the quark-gluon 
plasma, entropy production in the quark-gluon plasma, hadronization, and resonance 
production in the hadron phase and their decays. It is thus to be emphasized that 
the understanding of the final state observables requires the understanding of the 
whole period of the time evolution in heavy ion collisions. We will come back to this 
issue in Sec. [6j 

§3. Hydrodynamical expansion 

In the early stage of the RHIC operation, only hydrodynamical models were able 
to explain the strong elliptic flowp^ which was solid evidence of the creation of the 
strongly interacting QGP at RHIC. However, detailed analyses on experimental data 
eventually revealed the limitations of hydrodynamical models. They have difficulty 
in explaining, for example, the elliptic flow at forward/backward rapidity, HBT 
results, the centrality dependence of elliptic flow, and so forth. At the same time, 
viscosities became one of the most central topics in relativistic heavy ion collisions, 
while most of analyses had been carried out with ideal hydrodynamics. After the 
discovery of the strongly interacting QGP, the main interest at RHIC shifted to the 
understanding of detailed properties of the strongly interacting QGP. 

3.1. Hydrodynamical equations 

The basis of hydrodynamical models is the energy and momentum conservation, 

a At T^(x) = 0, (3-1) 

where T tiV {x) is the energy momentum tensor. In the case of ideal relativistic fluid, 
the energy momentum tensor is given by 

T» v (x) = [e(x) + p{x)]u^(x)u v (x) - p{x)g^ v , (3-2) 

where e(x), p(x), and u^(x) are the energy density, pressure, and four velocity, respec- 
tively. Eq. (|3-ip is solved numerically simultaneously with the charge conservation 
relation, 

d lx f(x) = 0. (3-3) 

When one starts to include the effects of dissipation into relativistic hydrody- 
namics, one is confronted with a rather complicated situation. One of the difficulties 
is that naive introduction of viscosities, first order theory (i.e., first order in gra- 
dients) suffers from acausality. The heat conduction equation allows instantaneous 
propagation of heat because of its parabolicity. The acausality of first order hy- 
drodynamics stems from the same reason. In order to avoid this problem, second 
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order terms in heat flow and viscosities have to be included in the expression for the 
entropy,^ '^'^'^'^'^'^D but the systematic treatment of these second order 
terms has not been established. Although there is remarkable progress toward the 
construction of a fully consistent relativistic viscous hydrodynamical theory, there 
are still ongoing discussions about the formulation of the equations of motion and 
about the numerical procedures.^ 1 

The basic tenet that has to be given up in dissipative hydrodynamics is the as- 
sumption of a uniquely defined local rest frame. Away from equilibrium the vectors 
defining the flows of energy, momentum, and conserved charges can be misaligned. 
We can still define a local rest frame by just choosing a velocity u^^x) in the labo- 
ratory frame. Then the energy-momentum tensor and the conserved charge current 
take more general involved forms, 

T fJ,u (x) = [e(x) + p(x) + II (x)]^ {x)u u (x) — \p(x) + nixW + 2W^u v) + 7T^ , 

(3-4) 

f(x) = n{x)u^ + V 1 , (3-5) 

where and are corrections to the flow of conserved charge and energy that 
are orthogonal to u^ 1 and T^ u u u , respectively, tt^ (with the orthogonality conditions 
v — 0) is the symmetric traceless shear stress tensor, and 77 is the 
bulk stress. (• • • ) indicates the symmetrization with regard to the indices. Usually 
is chosen to define one of the two standard frames: the Eckart frame where the 
velocity is given by the physical flow of net charge (then = 0), or the Landau 
frame where the velocity is given by the energy flow (then = 0). We refer the 
reader to the article by Muronga and Rischke for further discussions.^ 1 

At first order the new structures are proportional to gradients of the velocity 
field u^, and only three proportionality constants appear: the shear viscosity 77, the 
bulk viscosity £, and the heat conductivity k. With the usual definitions the first 
order relations in the Landau frame are^ 1 



n = -cv^ , (3-6) 

f? = -*^V"±, (3-7) 
e+p T 

t<v = 2? 1 V < ^u u> , (3-8) 

where = — (e + p)/nV >1 is the heat flow, V M = (g^ v — u fJ, u v )d u is the covari- 
ant derivative orthogonalized to the flow vector, T and [i are the temperature and 
chemical potential for the conserved charge, respectively. Here only one conserved 
charge is assumed, but the extension to cases with more than two conserved charges 
is straightforward. (• • • ) refers to the symmetrization of indices with the trace sub- 
tracted. The entropy current receives additional contributions beyond the equi- 
librium term su^ and one can show that all three transport coefficients are positive, 
demanding that the entropy is strictly non-decreasing, c^S^ > 0. 

At second order many more new parameters, related to relaxation phenomena, 
appear. Currently, most of viscous hydrodynamical calculations use the relativis- 
tic dissipative equations of motion that were derived phenomenologically by Israel 
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and StewarlP^ and variants of those, while some use the method by Ottinger and 
GrmelaP^'^'^ See for example Ref . ' '•' Recently, a second-order viscous hydrody- 
namics from AdS/CFT correspondence was derived,® as well as a set of generalized 
Israel-Stewart equations from kinetic theory via Grad's 14-momentum expansion, 
which have several new terms."*" On the other hand, however, a stable first-order 
relativistic dissipative hydrodynamical scheme was also proposed on the basis of 
renormalization-group consideration.®' 1 ^ 

In heavy ion physics, the shear viscosity, in particular its ratio with the entropy 
density, 77/s, has been attracted most of the attention among the above three trans- 
port constants. Interesting seminal investigations on the effects of bulk viscosity have 
begun,®'® while heat conductivity still has not been investigated systematically 
in connection with RHIC data. As for the second order parameters, there has been 
only little systematic investigations, either. 

Once equations of state are given, Eq. f)3- 1[) is solved, the effect of the phase 
transition being automatically taken into account, which is one of the advantages 
of the hydrodynamical model. We will discuss this feature in the next subsection. 
Recently the lattice (-inspired) equation of state which is connected to equation of 
state of resonance gas at low temperature is mostly used in hydrodynamical models. 
As the fireball expands, the temperature and density inside become so small that the 
assumption of the hydrodynamical picture becomes inapplicable any more. Thus, 
models with only a hydrodynamical component cannot describe the whole stages of 
the relativistic heavy ion collisions from the thermalized quark-gluon plasma stage to 
the kinetic freezeout. One way for such description is to connect a hydrodynamical 
simulation to a hadron based event generator as discussed in Sec. [6j 

3.2. Equation of State and Transport Coefficients 

One of the advantages of hydrodynamical models over other phenomenological 
models is their direct relation with the equation of state of QCD. Using the hydro- 
dynamical models one can find directly the consequence of the phase transition in 
experimental observables. In hydrodynamical models, historically, the equation of 
state with a first order phase transition based on the bag model widely has been 
used, because of its simplicity and lack of conclusive results on the equation of state 
of QCD. In recent hydrodynamical calculations, lattice(-inspired) equation of state 
has been widely employed, because of the development of thermodynamical analy- 
ses based on the first principle calculation, lattice QCD simulation. The equation of 
state of QCD for 2+1 flavors and also that including charm quark (2+1+1 flavors) by 
means of lattice simulations were reported by Borsanyi et al.® HoTQCD collabora- 
tion investigates chiral and deconfinement aspects of the phase transition with 2+1 
flavors using several kinds of staggered fermions.® There was some difference be- 
tween the two groups in the critical temperature, the temperature dependence of the 
so-called interaction measure and so on, but gradually the difference is disappearing. 
At the same time Wilson fermion simulation is also applied to the analysis of QCD 
thermodynamical properties.® Quantitative analyses of the QCD thermodynamics 
with the lattice simulation have just started. For conclusive results, improvement 
in actions, approach to continuum limit, simulation at the real pion mass and so on 
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need to be done. 

Simulations at finite chemical potential on the lattice suffer difficulty in execution 
of Monte Carlo simulations owing to the sign problem in the fermion determinant. 
Although there have been some developments recently such as the Taylor expan- 
sion method and reweighting method^ ^ ^ Mi Mi MM MM) Mb M the devel- 
opment in the study of finite chemical potential lattice simulations is relatively slow 
compared to lattice QCD study at the vanishing chemical potential. Though a lot of 
attempts to circumvent the problem have been made, we still need a breakthrough 
to explore the QCD properties in the whole region on the T-fi plane. 

In relativistic viscous hydrodynamical models, we need to input more informa- 
tion related to transport coefficients, shear and bulk viscosities, heat conductivity, 
and relaxation times, besides equations of state. The investigation of the transport 
coefficients of strongly interacting QGP and hadron gas is one of the most difficult 
problems in the field. There are several studies on transport coefficients with var- 
ious approaches; AdS/CFT,^ lattice QCD,^'^ finite temperature perturbative 
QCDpS>'E3'E3 microscopic transport models MjiEB ,EQ MbEB anc j re lativistic quan- 
tum Boltzmann approach. 841 However, to reach a conclusive result on the transport 
coefficients, in addition to theoretical calculations, it plays an important role to ex- 
tract transport coefficients from comparison of phenomenological model analyses and 
experimental data. 

3.3. Jet Energy Loss 

One of the most interesting features that are related to hydrodynamical expan- 
sion is the jet energy loss. A lot of interesting experimental results which suggest 
the existence of very large jet energy loss are reported. To explain these results 
interaction between jets and medium need to be understood. At the moment, there 
are perturbative QCD based approaches, the higher twist formalism, the AMY for- 
malism, the GLV formalism, and ASW formalism, and AdS/QCD approach for the 
jet energy loss (References are found, for example, in RefJ^SJ). 

Despite the large amount of effort put into the development of perturbative de- 
scription of the hadron production in heavy ion collisions, there are uncertainties 
remaining about the exact nature of jet-medium interactions in the kinematic and 
temperature regimes relevant at RHIC. As a whole, the above four approaches de- 
scribe RHIC data well, but they reach very different quantitative conclusions about 
the quenching strength, or transport coefficient, q. This does not come as a big 
surprise since the approaches differ in some of their basic assumptions, and there 
are large uncertainties in modeling hard probes beyond the calculation of the energy 
loss rate for a quark or gluon. 

Currently the big picture can be summarized as follows: perturbative calcula- 
tions under various assumptions are compatible with RHIC data, but the constraints 
are insufficient to rule out any of the models. The experimental constraints are 
also insufficient to completely exclude non-perturbative mechanisms for jet quench- 
ing. Calculations using the AdS/CFT correspondence to model strongly interacting 
QCEpJMJMI' can describe the same basic phenomenology. Most likely this challenge 
to perturbative QCD can only be answered at LHC. The extrapolation of jet quench- 
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Table III. pQCD-based energy loss models!^ This table summarizes some of the key assumptions 
of the four perturbative calculations discussed in the text. The models differ in the assumption 
on the medium (thermalized, perturbative), kinematics, scales (E = energy of the parton, kr = 
transverse momentum of the emitted gluon, /i = typical transverse momentum picked up from 
the medium, T = temperature, A — typical momentum scale of the (non-thermalized) medium, 
x = typical momentum fraction of the emitted gluon), and the treatment of the resummation. 



Model 


Assumption about the medium and kinematics 


Scales 


Resummation 


GLV 


static scattering centers (Yukawa), opacity 
expansion 


E > k T ~ n, x « 1 


Poisson 


ASW 


static scattering centers, multiple soft scat- 
tering (harmonic oscillator approximation) 


E > k T ~ n, x « 1 


Poisson 


HT 


observable matrix elements at scale A 
(thermalized or non-thermalized medium) 


E > k T > A ~ fi 


DGLAP 


AMY 


perturbative, thermal, g « 1 (asymptot- 
ically large T) 


£ > T > <?T ~ ft 


Fokker-Planck 



ing to larger jet energies is significantly different in strong coupling and perturbative 
scenarios.^ It is also possible to assume a small regime of strong non-perturbative 
quenching around T c together with the perturbative quenching at higher temper- 
atures. Such hybrid scenarios might be hard to distinguish experimentally. One 
of such pictures was recently explored by Liao and ShuryakP^ They found that 
a "shell" -like quenching profile in which quenching is enhanced around T c can give 
better simultaneous fits to single hadron suppression and elliptic flow. For more 
details, see, for example, RefP^ 

3.4. Hydrodynamical Expansion 

Before going to the discussion of experimental data on hydrodynamical expan- 
sion, in this subsection, we show the behavior of the temperature and chemical 
potential in hydrodynamical expansion. As an example, we depict the behavior of 
isentropic trajectories in the T — fi plane for Au+Au ^/snn = 200 GeV central col- 
lisions in Fig. HI The dotted line stands for the phase boundary between the QGP 
and the hadronic phase (Note that due to small baryochemical potentials, the phase 
boundary is an almost flat line at T c = 160 MeV). In addition to the central cell, 
we also investigate the isentropic trajectory of a cell close to the surface of the ini- 
tially produced QGP. Whereas the isentropic trajectory of the central cell located at 
(0, 0, 0) starts in the QGP phase (solid line), the cell at the initial surface of the QGP 
(dashed line) only exhibits an evolution from the mixed phase to the hadronic phase. 
Both trajectories are terminated at the freezeout temperature, Tf = 110 MeV. 
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Fig. 4. Isentropic trajectories on the 
T — (i plane in the case of the 1st or- 
der phase transition. Solid (dashed) 
line stands for the time evolution of 
the cell which is located at (a;, y, rj) — 
(0,0,0)(= (0,6,0)) at the initial time. 
The dotted line represents the phase 
boundary. 



10 20 30 
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3.5. Experimental Data and Discussion 

Calculated results shown in this subsection are based on RefP 1 To emphasize 
the importance of the final state interactions in the freezeout process in understand- 
ing hadron observables in relativistic heavy ion collisions, in this subsection we only 
show results of pure hydro dynamical calculation. The calculation is performed with 
a (3+1) dimensional ideal hydrodynamical model with the Glauber type initial con- 
dition and a bag model equation of state. Effects of final state interactions will be 
discussed in Sec. [H This calculation, thus, can be considered as a baseline for recent 
more realistic hydrodynamical models. 

In the initial energy density distribution, the maximum value of the energy den- 
sity is 55 GeV/fm 3 , which is relatively higher than in other hydrodynamical models, 
because in the pure hydrodynamical calculation we just use a single freezeout tem- 
perature and neglect both resonance decays and final state interactions. Usually 
parameters for the initial condition are set from comparison with experimental data 
of single particle distributions, rapidity distributions and Pt distributions in cen- 
tral collisions. Therefore hydrodynamical models have prediction power for other 
physical observables such as collective flow and the impact parameter dependence of 
various physical observables. We include the small baryon number density in these 
calculations. The starting time of hydrodynamical expansion is tq =0.6 fm. For 
the details of the parameters used in the calculations and the equation of state, see 
RefP 

First we show two examples which clearly show the limitation of pure hydrody- 
namical models. One is the Pt spectra of p and multistrange particles and the other 
is the elliptic flow as a function of the rapidity. Figure [5] shows the Pt spectra of 
7T, K, and p in Au + Au at y / SNN = 200 GeV for central collisions. Our calcula- 
tion succeeds in reproducing the ir spectra measured by PHENIXP^ up to Pt ~ 2 
GeV. However, due to the model assumption of chemical equilibrium down to the 
(low) kinetic freezeout temperature, we fail to obtain the correct normalization and 
hadron number ratios, even though the shape of the Pt spectra of p and multistrange 




Modeling a Realistic Dynamical Model 



13 



baryons (shown in Fig. [6]) is close to that of the experimental data. In order to obtain 
the proper normalization for the p spectra and hadron number ratios, we adopt a 
procedure outlined in Ref.p2i it renormalizes the proton Pt spectra using the p to 
7r ratio at the critical temperature. It is straightforward to extend this procedure 
to hyperons and multi-strange baryons as well, even though we choose to show the 
raw, unrenormalized, result for the multi-strange baryons in Fig. [6] to elucidate the 
situation prior to the renormalization. 

The need for renormalizing the p spectra suggests that the assumption of a per- 
sistent chemical equilibrium throughout the hadron phase until kinetic freezeout is 
not realistic and that an improved treatment of the freezeout process is required. One 
method to deal with the separation of chemical and thermal freezeout is the partial 
chemical equilibrium model (PCE)W^^^^^ below the chemical freezeout tempera- 
ture T c h one introduces a chemical potential for each hadron whose yield is supposed 
to be frozen out at that temperature. While the PCE approach can account for 
the proper normalization of the spectra, it fails to reproduce the transverse momen- 
tum spectra and mass dependence of the elliptic flow.^3 In Sec. [6j we shall utilize 
our hybrid hydro+micro model to decouple the chemical freezeout from the kinetic 
freezeout. In these hybrid approaches j®'^'® the freezeouts occur sequentially as a 
result of the microscopic evolution. Flavor degrees of freedom are treated explicitly 
through the hadronic cross sections in the microscopic transport, which leads to the 
proper normalization of all hadron spectra. 
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Fig. 5. Pt spectra for ir + , K + , and p in cen- 
tral collisions with PHENIX dataP^ For the 
proton we rescale our result using the ratio at 
the chemical freezeout temperature. See the 
text. 
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Fig. 6. Pt spectra for multistrange baryons in 
central collisions with STAR dataP^ In this 
pure hydro calculation no additional proce- 
dure for normalization has been done. 



Figure [7] shows the elliptic flow as a function of the pseudo-rapidity r) in central 
(3-15 %) and mid central collisions (15-25 %). In both cases our hydrodynamical 
model calculations overestimate the elliptic flow at forward and backward pseudo- 
rapidities, similarly to the results shown in RefP^ At large forward and backward 
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rapidities the assumptions of perfect hydrodynamical models such as local equilib- 
rium, vanishing mean free path, and negligible viscosity effect, are apparently no 
longer valid. The difference between experimental data and calculated results at 
forward and backward rapidities increases with the impact parameter, implying de- 
crease of the volume in which the hydrodynamical limit is achieved. 

Next we move on to the topic of jet energy loss in expanding hydrodynami- 
cal medium as one of interesting applications of hydrodynamical model. Figure [8] 
shows results from a comparative study by Bass et al. at RHICP® j e ts are propa- 
gated through a medium described by the same hydrodynamics, using three different 
schemes for energy loss: ASW, Higher Twist (HT), and AMY. In Fig. O-R^a as a 
function of Pt for two different centrality bins is shown. It shows that the in- 
dependence and the centrality dependence of Raa are described well by all three 
models. Each model has one free parameter that has been adjusted: the strong 
coupling as in AMY, % for the overall fit parameter in HT, K (q = K ■ 2 ■ e 3 / 4 ) 
in ASW. These parameters are fixed from the comparison with Raa data in central 
collisions (the top figure in Figj8|). 

This study confirms the remarkably large q in the ASW model compared to that 
in the HT approach. For the case where the quenching strength scales with e 3//4 the 
initial values q rRl found for the quark at the center of the fireball in a central collision 
are™ 

q = 18.5 GeV 2 /fm for ASW, q = 4.5 GeV 2 /fm for HT (3-9) 

and for the case where the quenching strength scales like the temperature T they 
are 

q = 10GeV 2 /fm for ASW, q = 2.3GeV 2 /fm for HT, q = 4.1GeV 2 /fm for AMY. 

(3T0) 

Note that the jet propagation in AMY model is calculated self-consistently as a 
function of the local temperature so that there is no difference between the two 
cases. 

This comparison is unique and very valuable in the respect that the same ini- 
tial hard cross sections and the same maps for the fireball, from the same (3+1)- 
dimensional ideal hydrodynamics were used. Any differences in the extracted values 
of q are solely due to differences in models for the jet energy loss. One of the conclu- 
sions is that our current knowledge applied to Raa leaves rather large uncertainty 
in the determination of q. 

Here we make a short comment on the jet energy loss at LHC. ALICE collabo- 
ration reported the nuclear modification factors Raa as a function of Pt of charged 
particles in central and peripheral Pb-Pb collisions at -^/saF/v = 2.76 TeVP^ They 
found that in central collision Raa is more suppressed up to Pt = 7 GeV compared 
to the PHENIX and STAR experiments at RHIC, which suggests that a much denser 
medium is produced and stronger parton energy loss occurs at LHC. Raa decreases 
with Pt for 2 < Pt < 7 GeV and it takes the minimum around Pt = 7 GeV and 

*' Note that the meaning of go in Eqs. (|3-9[l and (|3-10|l is different from the overall parameter 
go in HT. Both go are commonly used in the literature. 
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Fig. 7. Elliptic flow as a function of w with PHO- 
BOS experimental dateP* for central (3-15 %) 
and mid central (15-25%) collisions. Impact 
parameters are set to 4.5 (central) and 6.3 
(mid central) fm, respectively. 
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Fig. 8. Raa as a function of Pt for central (top) 
and mid-central (bottom) collisions calculated 
from ASW, Higher Twist, and AMY energy 
loss models. The single parameter in each 
model has been fitted to describe the data by 
PHENIXP* 



interestingly, increases with Pt for Pt > 7 GeV. The increase of Raa at high Pt 
was observed, for the first time, at LHC. It was not seen clearly at RHIC. 

§4. Numerical Schemes for Solving Hydrodynamical Equations 

In this section, we first explain the current status of hydrodynamical models 
from the point of view of numerical schemes for relativistic hydrodynamical equa- 
tions. As we showed in Tabs. U and HU hydrodynamical models are categorized 
into ideal and viscous hydrodynamical models. In addition, the difference of each 
hydrodynamical model is found in the space-time dimension of simulation, initial 
conditions, equations of state, and prescriptions for freezeout process. In current 
understanding, the most realistic hydrodynamical model should have the following 
features: viscosity effects, (3+1) dimensional space-time expansion, event-by-event 
fluctuated initial conditions, lattice QCD inspired equations of state, and freezeout 
process which is described by hadron based cascade models. On these issues ideal 
hydrodynamical models have been studied more deeply and its status is considered 
to be more mature than viscous hydrodynamical models. The investigation with 
viscous hydrodynamical models has just started. 

In addition to the above issues, an important ingredient in hydrodynamical mod- 
els should be taken into account seriously. It is what numerical scheme is adopted for 
solving the relativistic hydrodynamical equation. Up to now, only a little attention 
has been paid to this point. As long as we analyze multiplicities and collective flow 
using smoothed initial conditions, which numerical scheme to choose is not so im- 
portant. However, when we start to investigate viscosity effects and event-by-event 
fluctuations, we need to choose suitable numerical schemes carefully. 
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The relativistic hydro dynamical equation is a non-linear system equation, whose 
analytical solution is usually difficult to find. However, from the one dimensional 
wave equation dtu + cd x u = 0, which is much simpler than the actual relativis- 
tic hydrodynamical equation, we can explore a suitable numerical scheme by com- 
paring with the analytical solution. The naive differential scheme such as the 
FTCS (Forward-Time Central-Space) scheme = Uj — ^K n j+i — with 

v = causes an unphysical oscillation, which continues to grow after several time 
steps. In order to stabilize the unphysical oscillation in the FTCS scheme, one can 
use, for example, the Lax-Friedrich scheme. In this scheme on the right hand side 
of the FTCS scheme is replaced by the averaged value ^(u'j_ 1 + u"j +1 ). This scheme 
is stable if the Courant-Friedrichs-Lewy condition (CFL condition) [\u\ < 1) is sat- 
isfied, but it is known that it suffers huge numerical dissipation. In other words, 
the average manipulation introduces the artificial viscosity. One of the improved 
versions of Lax-Friedrich scheme is the Lax-Wendroff scheme, which has the second 
order accuracy in time and space, 



This scheme is stable, but the unphysical oscillation occurs at discontinuity. In order 
to avoid this problem introduction of artificial viscosity or flux limitation is needed. 
In conjunction with this fact we just cite Godunov's theorem: no second-order or 
higher order explicit monotonous scheme exists. Systematic discussion on numer- 
ical schemes is out of scope of this article. For details, please sec, for example, 
j^ e f[i02|i p rom these lessons for numerical schemes for the wave equation, we can 
deduce those for the relativistic hydrodynamical equation, though actual numerical 
tests are indispensable, i) First order accuracy scheme: For stability at discontinu- 
ity some average manipulation is needed, but this introduces large dissipation, ii) 
Second order accuracy: To remove numerical oscillation at discontinuity, one needs 
an artificial viscosity or a flux limiter. The former suggests that a simple central 
difference scheme of the first order accuracy might have a large artificial viscosity, 
which is crucial for the study of the viscosities of the matter created in relativistic 
heavy ion collisions. 

For event-by-event fluctuated initial conditions, shock- wave capturing schemes 
play an important role in dealing with discontinuity in the initial conditions. A lot 
of shock-wave capturing schemes have been proposed and developed. On the other 
hand, in relativistic heavy ion collisions, SHASTA, rHLLE, and KT algorithms are 
mainly used. In particular, SHASTA algorithm, which is widely used in the study 
of relativistic heavy ion collisions, is known as the first version of Flux Corrected 
Transport (FCT) algorithm In this algorithm, first a low-order solution is cal- 
culated. It incorporates large numerical diffusion effect. In the second step, as much 

*' The upper index (n) stands for the time step and the lower index (j) stands for the spacial 
position. 
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diffusion as possible is removed from the low-order solution. The amount of antid- 
iffusion fluxes is determined with the mask coefficient A^. This default value is 
^ad = 1) which can be set to lower values to reduce the amount of antidiffusion. 
In RefP® an interesting demonstration of the interplay between the numerical vis- 
cosity and physical viscosity was shown. A comparison was made with numerical 
solutions of the Riemann problem, one with the standard mask coefficient ^4 a d = 1.0 
SHASTA with a small physical viscosity r]/s = 0.01 and the other with a reduced 
mask coefficient = 0.8 SHASTA with vanishing viscosity. It was found that the 
difference of two numerical calculations is very small, which implies that it might be 
difficult to distinguish between the physical viscosity and the artificial dissipation. 
For quantitative investigation of the viscosities in the quark matter, one need to 
estimate the amount of artificial viscosity existing in numerical schemes carefully 
or choose numerical schemes known to have less artificial viscosity. Furthermore in 
Ref.p® a comparison among the different shock- wave capturing schemes, SHASTA, 
KT, and NT schemes was made and found that all the algorithms reproduce the 
analytic solution with nearly the same accuracy and numerical artifact (Figs. [9] and 

Ho]). 

Recently, RefU^ proposed a fast numerical scheme for causal relativistic hy- 
drodynamics with dissipation for analyses of relativistic high energy collisions, which 
is based on RefP^ In the numerical scheme, Israel-Stewart theory was adopted and 
the Israel-Stewart equations were decomposed into the in viscid part and the d i ssip a- 
tive part. For the inviscid part a relativistic Riemann so l ve ^MSMmBEmEmEmEmmBEm 
is used and for the dissipative part the Piecewise Exact Solution methocP^' is em- 
ployed in order to achieve less artificial dissipation and less computational time. 
Riemann solvers are numerical fluid dynamical solvers which are based on the con- 
cept of the Riemann problem. Several kinds of Riemann solvers have been proposed: 
Godunov method, Roe scheme, HLLE and HLLC solvers, and so on. Each solver has 
advantages and dis advantages in numerical cost, calculational accuracy, and coding 
complexity! 1 ^ 1 fiSD T n order to obtain the correct physical viscosity, for the inviscid 
part as well as the dissipative part, we need to choose numerical schemes which suf- 
fer as little artificial dissipative effect as possible. In Ref.p^ 1 for the inviscid part 
they use the relativistic Godunov method which is based on the exact solution of 
the Riemann problem. Figures l9l and [TUl show numerical solutions of the relativistic 
Riemann problem for ideal fluid on a grid with ^=100 cells with Ax=0.1 fm after 
Nt = 100 time steps at i = 4 fin together with the analytic solution. The initial 
conditions are set as the same as those in RefP® On the left the initial temperature 
is To = 0.4 GeV and on the right To = 0.2 GeV, and the energy density e is given by 
e = ^§T 4 (g = 16). The numerical results for SHASTA, KT, and NT schemes are 
taken from RefP^I SHASTA, KT, and NT schemes reproduce the analytic Riemann 
solution with nearly the same accuracy and numerical artifacts. On the other hand, 
we can see that the red lin c 1 105 ^ is closer to the analytical solution around x = 3 fm, 
which implies that this scheme suffers less artificial dissipative effect and is more 
suitable for physical viscosity analyses. A comparison of rHLLE and SHASTA was 
done in RefpM It was found that rHLLE has almost the same artificial viscosity as 
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SHASTA. 




Fig. 9. Energy density distributions as a function 
of x from numerical results (SHASTA (short- 
dashed line), KT (dotted line), NT (long- 
dashed line), and RefP 3 (red solid line)) 
together with the analytic Riemann solution. 



Fig. 10. Velocity distributions as a function of 
x from numerical results (SHASTA (short- 
dashed line), KT (dotted line), NT (long- 
dashed line), and RefpU (red solid line)) 
together with the analytic Riemann solution. 



§5. Hadronization 



In hydro dynamical models the hadronization process from the QGP phase to the 
hadron phase is naturally encoded in the equation of state. At RHIC hydrodynam- 
ical calculations do well for explaining a lot of experimental data, the Pt spectra, 
and the elliptic flow up to Pt ~ 2 GeV. However, above Pt ~ 2 GeV the results 
from hydro dynamical models start to show deviation from experimental data, which 
suggests that other mechanisms become dominant at higher Pt instead of the hy- 
drodynamical feature. In the intermediate Pt ( 2 < Pt < 6 GeV) the recombination 
mechanism is dominant and at high Pt ( > 6 GeV) the fragmentation mechanism 
plays an important role. These transverse momentum regions where the hydro- 
dynamical picture, recombination model, or fragmentation mechanism works well 
depend on collision energy. For example, at LHC the hydro dynamical description 
may work up to Pt ~ 4 — 5 GeV.^ There have been attempts to use thermodynam- 
ical properties which are described with a hydrodynamical model as inputs of the 
recombination model and fragmentation mechanisml^^'G^ \i may become possible 
to construct more realistic dynamical models to understand the physics of relativis- 
tic heavy ion collisions by employing the recombination model and fragmentation 
process for hadronization in a hydrodynamical model. In the next subsection, we 
give a brief description of the recombination model. 

5.1. Recombination Model 

Quark recombination or coalescence is the best candidate for the physical mech- 
anism to explain a large amount of experimental data at the intermediate Pt at 
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RHIC. Now not only at LHC but also at relatively low collisions energy experi- 
ments, analyses with the recombination picture are ongoing. 121 ' First the quark 
number scaling law was found in the behavior of elliptic flow as a function of the 
transverse momentum, which is considered as a signature of the quark recombina- 
tion. For a quark phase with the elliptic flow u|(Pt) a t the time of the hadronization 
simple instantaneous recombination models predict 

v H(P T ) = nvl[^\ , (5-1) 

where n is the number of the valence quarks. This scaling law describes the key 
feature of experimental data at the intermediate Pt notably well. 

Generally, recombination models assume a universal phase space distribution of 
quarks at hadronization. Quarks turn into baryons, qqq — > B, and mesons, qq — > M. 
These processes are descr ibed either by using instantaneous projections of quark 
states onto hadron st at es >EH ,li25p ,|126J| or dy nam ical coalescence processes 
with finite width hadrons governed by rate equations Note that usually only 
the valence quarks of the hadron are taken into account although generalization has 
been attempted P^ 

The original instantaneous projection models explicitly preserve only three com- 
ponents of the energy-momentum four-vector in the 2 — > 1 and 3 — > 1 processes. The 
yield of mesons can be expressed through the convolution of the Wigner function 
W a t, for a pair of partons, a and b and the Wigner function <Pm encoding the meson 
wave function 

dN M f d 3 R fd 3 qd 3 r Txr r P „rP \ x , . 

-#p=] W?^j I'Y + q 2'2-- q ) $M(r ' q) ' 

(5-2) 

The quark Wigner functions are usually approximated by the classical phase space 
distributions. Hadron spectra at intermediate Pt are described well by considering 
the factorization of W a & into the thermal quark distributions 

W a6 (ri ! pi;r 2 ,p 2 ) = /o(ri,pi)/b(r 2 ,p 2 ) • (5-3) 

Correlations between quarks can be introduced to model correlations found between 
hadrons^ 2 ^ without interfering with the excellent description of the spectra and 
hadron ratios. 

The strength of the quark recombination picture is in its predictive power, which 
originates from its explaining all measured hadron spectra at the intermediate Pt 
basically with one parameter for the quark distribution function at hadronization. 
It has been shown that at low momenta resonance recombination is compatible 
with hydrodynamics and kinetic equilibrium but on the other hand, because 
thermalized states do not retain memories of previous time, all phenomena with 
long time scale in the equilibrated region should be expected to be described by 
hydrodynamics. They include the quark number and kinetic energy scaling observed 
at RHIC at low moment a P^ The possibility of including quark recombination 
explicitly into hydrodynamical model has been studied in RefsP^^^ 



20 



C. Nonaka and M. Asakawa 



5.2. Experimental Data and Discussion 

Using the recombination model, we can explore the production of quark soup, 
which is a consequence of the QCD phase transition. After the success of the recom- 
bination model at RHIC, the scaling property of the elliptic flow has been tested in 
a wide range of collision energy to investigate the properties of the strongly inter- 
acting QGP. The elliptic flows of ^,p, and A were measured in Pb+Pb collisions at 
^/sNN : = 17 GeV at SPS by NA49.^ They found that the quark number scaling 
in the elliptic flows of these particles holds in the Pt range covered by the data (up 
to Pr/n ~ 1 GeV. Fig. 6 in RefpSHJ). Since the Pt is limited to rather low values, 
however, this should not be viewed as a conclusive test for the quark recombination 
mechanism at SPS. 

At LHC the analyses of the quark number scaling in the elliptic and triangular 
flows of identified particles have just started. In Ref.p^ 1 ALICE shows the elliptic 
and triangular flows per constituent quark of tt^ , , and p as a function of the 
transverse kinetic energy (KEt) per quark for more central (10 -20 %) and more pe- 
ripheral (40-50 %) Pb+Pb collisions at ^/snn = 2.76 TeV. They report that within 
errors the flows of ir and K follow the scaling, while the flow of p deviates in the more 
central and more peripheral events. In the triangular flow they find the same feature 
as the elliptic flow, i.e., the triangular flow of p shows deviations from the KEt 
scaling. However, this investigation is done in a relatively low transverse kinetic 
energy region (KEt/u < 1 GeV) where the hydrodynamical behavior is expected 
to be dominant. Therefore, this deviation from the KEt scaling in the elliptic and 
triangular flows could be explained by the mass splitting and mass ordering real- 
ized by hydrodynamical motion. To obtain a conclusive result for the quark number 
scaling in the collective flow at LHC measurement at higher transverse kinetic en- 
ergy, where the effects of the recombination mechanism are more clearly observed, 
is indispensable. 

For the experimental confirmation of the quark number scaling in the elliptic 
flow, no particles are more important than the (j> meson. Because it is composed 
of one quark and one antiquark and its mass is close to that of the proton, anal- 
yses of the elliptic flow of <j> reveal the dominant process in the hadronization at 
the intermediate Pt'- recombination mechanism or thermal process as described by 
hydrodynamicsP^ZJ i n addition, because multi-strange hadrons have large mass and 
small hadronic cross sections, they should be less sensitive to hadronic rescattering 
in the later stage of collisions and a good probe of the early state of the collisions 
such as the partonic elliptic flow. STAR Collaboration measured the V2 of multi- 
strange hadrons (</>, E, and Q) in Au+Au collisions at ^Jsnn = 200 GeV and showed 
the quark number scaling works in the elliptic flow of the multi-strange hadrons.' 1 ^ 3 -' 1 
This implies that the partonic collectivity is built up at the top RHIC energy. 

Furthermore STAR Collaboration investigates the energy dependence of the 
elliptic flow as a function of the transverse momentum for <fi at RHIC. 134 ' The 
quark number scaling in the elliptic flows of ir, K, p, A, and E works well at 
both t/saF/v = 39 and 11.5 GeV in Au+Au collisions. They found that while at 
\J s nn = 39 GeV the elliptic flow of <j> also follows the quark number scaling, at 
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y / s]v r /v = 11-5 GeV the elliptic flow per constituent quark of <j> is much smaller than 
those of the other particles. Only the <fi meson is out of the quark number scaling 
of the elliptic flow at this energy. They conclude that it indicates the dominance 
of hadronic interactions with the decrease of the beam energy. To confirm this, 
measurement of the elliptic flow of Q would be helpful. 

At the top RHIC energy, the elliptic flow of hadrons which contain only light 
flavor quarks (u, d, and s) follows the quark number scaling. Does the quark number 
scaling work in the elliptic flow of hadrons which are composed of heavy flavors? Re- 
cent STAR measurement of the elliptic flow of J ftp gives the answer to this question. 
It shows that the elliptic flow of J ftp is consistent with zero over the measured Pt 
range, which is significantly smaller than that of cj) and inclusive charged hadronsP^ 
The picture that J /tp production is dominated by charm quark recombination with 
significant charm quark flow is disfavored. On the other hand, the elliptic flow of 
open charm mesons is well understood by the recombination of a light (anti) quark 
and an (anti)charm quark with elliptic flowP^SJ 

Next we argue that the recombination model is one of powerful tools to under- 
stand hadron properties and the QCD phase diagram with high energy heavy ion 
collisions. The investigation of the elliptic flow of hadron resonances with the recom- 
bination picture at RHIC makes it possible to know final state interactions in the 
freezeout processes and the structure of resonances In principle, there are two 
different mechanisms that contribute to resonance production and as a result there 
are two types of resonances: (1) primordial resonances - resonances produced from 
hadronizing quark gluon plasma (QGP mechanism), and (2) secondary resonances - 
resonances produced in the hadronic final state via hadron-hadron scattering (HG 
mechanism). For example, in the case of K$, which is composed of the d and s 
quarks, Kq are produced by the recombination of the d and s quarks in the QGP 
mechanism, while it is created via scattering of K and tt in the HG mechanism. 
The recombination model tells us that the scaling constant n of the elliptic flow of 
a hadron species is given by the number of constituent quarks which participate in 
the production of the hadron: in the case of the QGP mechanism n = 2, and in the 
case of HG mechanism n = 4. From these two contributions the measured elliptic 
flow is given by, 

^measured = + (1 _ r (p r ))^G > (5.4) 

where r{Px) is the fraction of resonances which were produced at the hadronization 
and whose decay products escaped from the hadron phase without rescattering. 
With t{Pt) one can investigate the cross sections of resonances in the hadronic 
medium. STAR attempted to extract the scaling constants (or apparent constituent 
quark number) from the measured elliptic flow of K°* to investigate its production 
mechanism. Unfortunately, the obtained value of n is n = 3 ± 2p^ from which we 
cannot draw a definite conclusion. This method is easily extended to explore the 
structure of exotic hadrons, such as tetraquark and pentaquark. Utilizing the scaling 
law of the elliptic flow, one can obtain the information on the structure of exotic 
hadrons, whether they have compact structure or are more like molecular bound 
statesPB 
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We can explore consequences of diquark and quark-antiquark clustering above 
the deconfinement temperature with event-by-event net charge fluctuations!^® Re- 
cently, lattice QCD calculations show that charmonia survive even well above the 
critical temper at ure >UM$ UMi w hich suggests the possibility that hadroniza- 
tion occurs via recombination of not only single q or q but also qq or qq. The D 
measure, the net charge fluctuation normalized by the entropy is defined by 

D = A((AQ) 2 )/N ch , (5-5) 

where ((AQ) 2 ) denotes the event-by-event net charge fluctuation within a given 
rapidity window Ay and iV c h is the total number of charged particles emitted in this 
window. For a free plasma of quarks and gluons D ~ 1, whereas for a free pion gas 
D ~ 4. However, experimental data at RHIC is rather close to the value of hadron 
gas D = 2.8 ± 0.05 in central Au+Au collisions at ^/s]^ = 130 GeV (STAR^ 
and D « 3 (PHENIX) i n the recombination scenario, the fluctuation of the net 
charge Q = YliQi n i 1S given as 

W 2 > = (Q 2 ) - (Q) 2 = £>) 2 (n*> + Y. c ?k(^)( n k)Mk, (5-6) 

i i,k 

(2) 

where cV are the normalized two-particle correlation functions. In the absence of 
two-particle correlations, Eq. (|5-6j) can be rewritten as (5Q 2 ) = l(N u + N u ) + l{N d + 
Nj+Ns+Ng), where N{ = (rii) denotes the average number of the constituent quarks 
with flavor i. Together with the statistical hadronization model and the multiplicity 
of experimental value, the net charge fluctuation in the quark recombination scenario 
is d(5Q 2 ) q /dy = 331 ± 27, which is close to the experimental value d(5Q 2 )h ax i/dy = 
368 ± 33p^ Furthermore the difference of the net charge fluctuations between the 
quark recombination and the experimental data suggests necessity of improvement 
on both theory and experiment sides. On e of such possibilities is to include qq and 
qq pairs to the hadronization mechanism J^SJ If the qq pair and qq are taken into 
account, Eq. (|5-6p is extended to 

W 2 > = ^2(Qi) 2 (Ni + Nj) + Y,(Qi + qj) 2 (nij) + ^2(qi ~ q 3 ) 2 (nij) , (5-7) 

i ij ij 

where the average numbers of diquarks and qq are proportional to the products of 
the individual quark numbers: (riij) = a(NiNj + NjNj), (fiij) = (3NiNj with relative 
pairing weights a and f3. They showed that experimental data can be fitted with an 
appropriate cho ice of a and (3. 

In Ref.p^J the problem of the gluonic degrees of freedom was also discussed. 
According to recent lattice calculations pSI4SDJfi2I the phase transition at small baryon 
chemical potential is crossover. This means that the hadronization at RHIC and LHC 
does not take place from a system where quarks and gluons are active degrees of 
freedom, but from a state with reduced entropy density. At the moment, no method 
on the lattice has been found to figure out what is the active degrees of freedom in the 
crossover region. In RefP^t it was shown that entropy density data on the lattice 
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are not inconsistent with that of a gas of constituent quarks, which is assumed in 
the recombination model. The success of the recombination model strongly suggests 
that the gluonic degrees of freedom disappear first when the temperature is decreased 
from well-above the (pseudo) critical temperature. However, it is needless to say that 
it is important to find a way to figure out how the active degrees of freedom change 
across the crossover region on the lattice in order to steady this picture. 

§6. The Freezeout Process 

As fireballs expand, the temperature and density inside them become small. Fi- 
nally the mean free path among particles inside the fireballs becomes so large that 
the assumption for the hydrodynamical picture becomes invalid. Currently two sep- 
arate freezeout processes are believed to occur successively in heavy ion collisions. 
One is chemical freezeout, at which the ratios of hadrons are fixed, and the other 
is thermal (kinetic) freezeout, at which the particles stop interacting. Recently final 
state interactions between the two freezeout processes are also included in hydrody- 
namical models by connecting hadron based event generators to the hydrodynamical 
calculations. In this picture, kinetic freezeout is found to be not an instantaneous 
process but a continuous one as we show in the following. 

6.1. Chemical Freezeout and Thermal Freezeout 

The chemical freezeout temperature and baryon chemical potential are extracted 
with the statistical model on the basis of the grand canonical formalism. Sur- 
prisingly, statistical models are in excellent agreement with experimental data for 
hadron ratios in a wide range of collision energy from SIS to RHICP^ At present 
the (pseudo) critical temperature suggested from the latest lattice calculation is 
T c = 150 ~ 160 MeV (see Subsec. I3.2[) . which is relatively low from the point of 
view of the statistical modelP^' For example, using a statistical model, Andronic 
et al. extract the chemical freezeout temperature T c ^ = 160 ~ 166 MeV for LHC 
Pb+Pb collisions at - v /snn = 2.76 GeV. It might be too early to take lattice results 
at the face value, but this suggests that the physical meaning of the temperature 
and chemical potential obtained by the statistical model need to be reexamined. At 
the same time, one needs to make it clear why the statistical model works very well 
not only in a wide range of collision energies from SIS to LHC but also for smaller 
systems such as p+p collisions 

At the thermal freezeout temperature, the mean free paths of the hadrons have 
grown to the order of the system size and hydrodynamical description, which requires 
very small mean-free path, is clearly no longer applicable. A first naive guess for 
the kinetic freezeout temperature Tf would be of the order of the pion mass (~ 140 
MeV) . Practically the value of the kinetic freezeout temperature in hydrodynamical 
models is determined from comparison with data of the Py-spectra. 

The task at the end of a hydrodynamical calculation is to populate fluid cells 
with particles with the final temperature and flow. For calculations of single particle 
spectra, the simple assumption of a sudden freezeout process at a certain proper time 
for each fluid cell is often adopted, neglecting the reverse process from particles to 
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the hydro dynamical medium. Under this assumption the Cooper-Frye formular^ is 
widely used, 

e #p =T,jBrL d(TMjPM exp[(p^- / , / )/r f ]±i' (6 ' 1} 

where is the degeneracy factor of hadron h and Tf and //f are the freezeout temper- 
ature and baryon chemical potential, respectively, and +(— ) is for fermions (bosons). 
d(Tu is obtained by calculating the normal vector on the freezeout hypersurface E. 
Once these quantities are given, using Eq. (|6-1|) . one can calculate the distribution 
of any particle after freezeout. 

More realistic models have been investigated. One of them is the Continuous 
Emission Model (CEM), in which particles are emitted continuously from the whole 
expanding volume of the system at different temperatures and different timesP^ 
In the early days of hydrodynamics only kinetic freezeout was implemented. In- 
deed, at lower collision energies such as at SIS and AGS, the separation between 
the chemical freezeout and kinetic freezeout points is not large on the T — /i plane. 
However, at RHIC there is a significant difference between kinetic freezeout temper- 
atures from hydro-inspired models and the chemical freezeout temperature from the 
statistical modelP^ This fact also manifests itself through the failure to get the 
correct absolute normalization of some P^-spectra, e.g., that of protons in hydrody- 
namical calculations, with only a kinetic freezeout For these reasons, attempts 
to model separate freezeout processes consistently via modified equations of state 
were started P'^^'EU 

6.2. Final State Interactions 

It turns out that some experimental data are still not understood in a satisfactory 
way even with the two separate freezeouts. For example, mean transverse momentum 
(Pt) as a function of particle mass does deviate from the linear scaling law, which 
suggests significant final state interactions in the hadronic phase. 81 To explain these 
effects, and to account for the apparently large viscosities in the hadronic phase, 
as discussed before, hydro+cascade hybrid models were introduced. They use a 
hydrodynamical computation of the expansion and cooling of hot QCD bulk matter, 
and then couple the output consistently to a hadron-based transport model in order 
to take account of the final state interactions. A pioneering work on hydro+cascade 
hybrid models was done by Bass et alPSSl using UrQMD. Similar investigations 
were carried out in RefsP)® Fi gure [11] shows a schematic sketch of the hydro + 
UrQMD model. 8 ' At the switching temperature Tgw, the mean free path of hadrons 
becomes so large that the hydrodynamical picture becomes inapplicable. Thus after 
this point, the motion of hadrons is described by UrQMD. In practice, at Tsw hadron 
distributions are calculated with the Cooper-Frye formula in the fluid and the initial 
state of UrQMD is produced with them through the Monte Carlo methods. Then the 
UrQMD simulation is started. The reverse process from UrQMD to fluid dynamics 
is neglected in the simulation. 
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Fig. 11. Schematic sketch of the 3D hydro+UrQMD model. T c and and Tsw are taken to be 160 
MeV and 150 MeV, respectively. 

6.3. Experimental Data and Discussion 

The main purpose of our hydrodynamical model + hadron based event generator 
is to handle realistic final state interactions in the freezeout process. We need to 
choose an appropriate hadron based event generator for it. Hadron based event 
generators have been developed to understand experimental data especially at lower 
energy heavy ion collisions. Here we use the UrQMD model for the description of 
final state interactions in the hadron phase. It gives us reasonable hadron yields, 
single particle spectra and so on, and it is often used for understanding the baseline 
of experimental data of relativistic heavy ion collisions. 

Here we make a comment on the applicability of cascade models for the de- 
scription of the afterburner of hydrodynamical models. Recent lattice calculations 
(RefP^ 1 ) suggest that the hadron phase is well-described as a hadron resonance gas 
up to the vicinity of the phase transition. Thus, cascade models are, in general, ideal 
machinery to investigate its time evolution and freezeout process. 

If the final state interactions are implemented properly, hydrodynamical mod- 
els acquire more predictive power for experimental observables. We will show this 
by comparing single particle spectra and elliptic flows with experimental data. In 
this subsection we show results calculated with the same initial conditions and the 
equation of state as those in pure hydrodynamical calculations: Glauber type initial 
conditions and a bag model equation of state. In this case the maximum value of 
initial energy density is 40 GeV/fm 3 , which is smaller than that of the pure hydrody- 
namical model, because the inclusion of resonance decays and final state interactions 
to hydrodynamical models increases particle multiplicity. 

First let us consider the argument on the multiplicity, which we mentioned in 
Sec. We take tt + as an example, because it is one of the most dominant particles 
in the charged hadron multiplicity and affected a lot from resonance decays. Figure 
[T2l shows the pseudo rapidity distributions of n + from the fluid at T sw (solid line) , 
that from the fluid plus resonance decay effects (open circles), and its final state 
multiplicity (solid circles). We find that the resonance decay effect is huge and that 
the multiplicity including resonance decay effect (open circles) becomes almost twice 
as many as that of hydro at T sw (solid line). Furthermore, the additional gain from 
final state interactions is seen. This result clearly shows that the dynamics of the 
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hadron phase such as resonance decays and final state interactions, is important also 
in high energy heavy ion collisions. In short, final state multiplicities cannot be 
predicted by determining only the initial state. 

Figure [131 shows the Pt spectra of ir + , K, and p in central collisions at ^/snn = 
200 GeV. The most remarkable feature, compared to the pure hydro dynamical cal- 
culation, is that the hydro+micro approach is capable of accounting for the proper 
normalization of the spectra for all hadron species without any additional correction 
as is performed in the pure hydrodynamical model. The introduction of realistic 
freezeout process provides therefore a natural solution to the problem of the sep- 
aration of the chemical and kinetic freezeouts in pure hydrodynamical calculation. 
Similar results have been obtained previously in 1+1 and 2+1 dimensional imple- 
mentations of the hydro+micro approach >II£6I 
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In Fig. [T3] we analyze the Pt spectra of multi-strange particles. Our results 
show good agreement with experimental data for A, S, and Q for centralities 0-5 
% and 10-20 %. In this calculation no additional procedure for normalization is 
needed. Recent experimental results suggest that at the thermal freezeout multi- 
strange baryons exhibit less transverse flows and higher temperatures closer to the 
chemical freezeout temperature compared to non- or single-strange baryons 
This behavior can be understood in terms of the flavor dependence of the hadronic 
cross section, which decreases with increasing strangeness content of the hadron. 
The reduced cross section of multi-strange baryons leads to a decoupling from the 
hadronic medium at an earlier stage of the reaction, allowing them to provide in- 
formation on the properties of the hadronizing QGP less distorted by hadronic final 
state interactions It should be noted that the analogous behavior has 
already been observed in experiments at CERN SPSP^O 
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Fig. 14. Calculated Pt spectra of multi-strange particles at centralities 0-5 % and 10-20 % and 
STAR dataP 3 



In Fig. [TJ] the mean transverse momentum (Pt) as a function of the hadron 
mass is shown. Open symbols denote the value at T sw = 150 MeV, corrected for 
hadronic decays. Not surprisingly, in this case the (Pt) follows a straight line, sug- 
gesting a hydro dynamical expansion. However, if hadronic rescattering is taken into 
account (solid circles), the (Pt) does not follow the straight line any more. The 
(Pt) of pions is actually reduced by hadronic rescattering (they act as a heat-bath 
in the collective expansion), whereas protons actually pick up additional transverse 
momentum in the hadronic phase. RHIC data by STAR collaboration are shown by 
the solid triangles - overall, the proper treatment of hadronic final state interactions 
significantly improves the agreement of the model calculation with the data. We 
note that our results confirm those previously obtained in 1+1 and 2+1 dimensional 
implementations of the hydro+micro approach]^ 'L^* demonstrating the robustness 
of the hydro+micro approach across the three different implementations of the hy- 
drodynamical and macro to micro transition components of the model. 

Let us now investigate the effect of resonance decays and hadronic rescattering 
on the pion and baryon transverse momentum spectra. Figure [16] shows the Pt 
spectrum of ir + at T sw = 150 MeV (solid line, uncorrected for resonance decays) 
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as well as the final spectrum after hadronic rescattering and resonance decays, la- 
beled as Hydro+UrQMD (solid symbols). In addition, the open symbols denote a 
calculation with the resonance decay correction performed at T sw , which we label 
as Hydro+decay. The difference between the solid line and open symbols therefore 
directly quantifies the effect of resonance decays on the spectrum, which is most 
dominant in the low transverse momentum region Pt < 1 GeV. Furthermore, the 
comparison between the open symbols and solid symbols quantifies the effect of 
hadronic rescattering: pions with Pt > 1 GeV lose momenta via these final state 
interactions, resulting in the steeper slope. 
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Figure [T7] shows the elliptic flow as a function of rj: the pure hydrodynamical 
calculation is shown by the solid curve, the hydrodynamical contribution at T sw is 
denoted by the dashed line and the full hydro+micro calculation is given by the solid 
circles, together with PHOBOS data (solid triangles). The shape of the elliptic flow 
in the pure hydrodynamical calculation at T sw is quite different from that of the 
full hydrodynamical one terminated at the kinetic freezeout temperature, 110 MeV. 
Apparently the slight bumps at forward and backward rapidities observed in the 
full hydrodynamical calculation develop in the later hadronic phase, since it is not 
observed in the calculation terminated at T sw . Evolving the hadronic phase in the 
hydro+micro approach increases the elliptic flow at the central rapidities, but not 
in the projectile and target rapidity regions. As a result, the result for the elliptic 
flow in the hydro+micro approach is closer to the experimental data, compared to 
the pure hydrodynamical calculation. 

In Fig.[18]we investigate the effect of hadronic rescattering on the duration of the 
freezeout process by comparing the calculation terminated at T sw without hadronic 
rescattering (open symbols) and the one including the full hadronic final state in- 
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teractions (solid symbols). If we terminate time evolution at T sw , most hadrons 
freezeout around tj = 10 fm/c, reflecting the lifetime of the deconfined phase (the 
tails of the distribution come from the decays of long-lived resonances). The inclu- 
sion of hadronic rescattering shifts the peak of the freezeout distribution to larger 
freezeout times (77 ~ 20 — 30 fm) , which provides an estimate on the lifetime of the 
hadronic phase around 10-20 fm. 
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Fig. 18. Freezeout time distribution of mesons 
for hydro+decay (open symbols) and hydro + 
UrQMD (solid symbols) at the mid-rapidity in 
the case of central collisions. 



The findings discussed in the context of the previous figure for pions and kaons 
are confirmed by analyzing baryons in the same way, which is shown in Fig. [19l here 
the top frame contains the analysis terminated at T sw and the bottom frame contains 
the calculation including full hadronic rescattering. Figures [TH] and [T2] suggest that 
the assumption of sudden freezeout, which is often used, for example, in the blast 
wave model, is not realistic under the existence of final state interactions. 



§7. Summary 



In this article, we outlined key issues in constructing a realistic and compre- 
hensive dynamical model for the description of the whole stages of relativistic high 
energy heavy ion collisions: initial conditions, hydrodynamical expansion, hadroniza- 
tion, and freezeout processes. Hydrodynamical models are a promising starting point 
for building multi-module models. One of the reasons is that they give global un- 
derstanding of experimental data at RHIC and their success at RHIC makes us to 
expect that hydrodynamical analyses will serve as the baseline for investigation also 
at LHC. In addition, one can easily input up-to-date knowledge on each stage in 
the time evolution to hydrodynamical models easily: more realistic initial conditions 
which contain even- by-event fluctuations, latest equations of state and transport co- 
efficients from lattice QCD, realistic hadronization mechanisms such as those from 
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the recombination model and fragmentation mechanism, and realistic freezeout pro- 
cesses including final state interactions. 

Presently, both ideal and viscous hydrodynamical models are utilized for the 
investigation of data obtained at RHIC and LHC. In the multi-module modeling, 
however, ideal hydrodynamical models have been studied more deeply and its status 
is considered to be more mature than viscous hydrodynamical models. However, 
gradually the majority of hydrodynamical models will become viscous hydrodynam- 
ical models instead of ideal hydrodynamical models, because the study of viscos- 
ity effects in physical observables and the estimate of the viscosities in QGP from 
hadronic observables are among the hottest topics in the field; inclusion of viscosities 
in hydrodynamics is inevitable for such investigation. At the moment, investigation 
with viscous hydrodynamical models has just started. Whereas their development 
is fast and remarkable, a lot of issues remain to be considered, implemented, tested, 
and improved. 

One of them is numerical schemes for solving the relativistic hydrodynamical 
equation, to which only a little attention has been paid up to now. In particular, the 
estimate of numerical viscosity is crucial for the study of the viscosities of the matter 



Modeling a Realistic Dynamical Model 



31 



created in relativistic heavy ion collisions. In Sec. 01 we showed that SHASTA, KT, 
and rHLLE schemes, which are mainly used in studies of high energy heavy ion 
collisions, have nearly the same accuracy and numerical artifact. Furthermore, we 
showed the results of shock tube test with a new numerical scheme which suffers less 
artificial dissipative effect and showed that it is more suitable for analyses of physical 
viscosities than SHASTA, KT, and rHLLE schemes. 

In summary, in this article we discussed essential ingredients in understanding 
entire stages of high energy heavy ion collisions in detail, together with interpretation 
of experimental data at SPS, RHIC, and LHC. We hope that the road map shown 
in this article will serve as a guideline for modeling a realistic dynamical model for 
high energy heavy ion collisions. 
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